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Abstract —Contemporary electricity distribution systems are 
being challenged by the variability of renewable energy sources. 
Slow response times and long energy management periods 
cannot efficiently integrate intermittent renewable generation 
and demand. Yet stochasticity can be judiciously coupled with 
system flexibilities to enhance grid operation efficiency. Voltage 
magnitudes for instance can transiently exceed regulation lim¬ 
its, while smart inverters can be overloaded over short time 
intervals. To implement such a mode of operation, an ergodic 
energy management framework is developed here. Considering 
a distribution grid with distributed energy sources and a feed- 
in tariff program, active power curtailment and reactive power 
compensation are formulated as a stochastic optimization prob¬ 
lem. Tighter operational constraints are enforced in an average 
sense, while looser margins are enforced to be satisfied at all 
times. Stochastic dual subgradient solvers are developed based 
on exact and approximate grid models of varying complexity. 
Numerical tests on a real-world 56-bus distribution grid and the 
IEEE 123-bus test feeder relying on both grid models corroborate 
the advantages of the novel schemes over their deterministic 
alternatives. 

Index Terms —Energy management, reactive power compen¬ 
sation, active power curtailment, stochastic optimization, dual 
decomposition. 


I. Introduction 

Distributed generation and the prospective integration of 
electric vehicles and elastic loads create unseen operational 
scenarios for distribution grids 0. Several utilities in the US 
currently experience issues with integrating residential- and 
commercial-scale solar generation. For example, solar farms 
oftentimes connected at the end of a long distribution feeder in 
distant rural areas, are routinely reported to introduce voltage 
regulation problems to the residential buses across the feeder. 
These frequently reversing power flows strain the apparent 
power capabilities of substation transformers. Moreover, data 
collected from residential PVs reveal that solar generation 
can fluctuate by up to 15% of their nameplate ratings within 
one-minute intervals m. The aforementioned issues critically 
challenge energy management of distribution grids. 
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On the other hand, contemporary distributed generation 
units are equipped with the so-termed smart power inverters 
that have two-way communication and computing capabilities, 
and thus offer unprecedented control opportunities 0. Lever¬ 
aging smart inverters for joint reactive power compensation 
and active power curtailment to achieve various objectives 
(power loss minimization, conservation voltage reduction, or 
voltage regulation) is considered here. Traditionally, distribu¬ 
tion grid voltage regulation is performed via load-tap-changing 
transformers, capacitor banks, and voltage regulators S). This 
utility-owned equipment involves discrete control actions, and 
its lifespan is affected by frequent switching operations 0, 
0- Regulating voltage under increasing generation from dis¬ 
tributed renewable sources would require even more frequent 
switching actions and perhaps further installations. 

In this context, recent research efforts have focused on 
engaging smart inverters in the energy management system 
(EMS) of distribution grids 11 j, 0, J7); especially, given 
that these inverters come with PVs and electric vehicles 
anyways. Customer-owned power inverters can be commanded 
to adjust reactive power injections within milliseconds and in 
a continuously-valued manner 0, 0. Albeit currently pro¬ 
hibited by some standards (see e.g., j9j), controlling reactive 
power through smart inverters has been reported to improve 
grid’s voltage profile, or even displace utility-owned voltage 
regulating equipment at more than 50% solar penetration 0. 

Using approximate grid models, voltage regulation is ef¬ 
fected through a multi-agent scheme in fTOl . and local control 
algorithms are devised in 0. Building on the exact full 
AC grid model, reactive power control is an instance of the 
optimal power flow (OPF) problem, which is non-convex in 
general ED- Recently, different convex relaxations have been 
proposed; see m for a review. In radial networks, the OPF 
can be relaxed into a second-order cone program (SOCP) 
via either polar coordinates ED, or the branch flow model 
D3; or into a semidefinite program (SDP) m. Although 
the two relaxations have been shown to be equivalent, M 
advocates using the SOCP one due to its lower computational 
complexity. Leveraging the SOCP relaxation, a two timescale 
conventional and inverter-based reactive power control has 
been formulated in na. Accounting for stochasticity, an 
adaptive local control algorithm for single-branch grids is 
developed in lfl8l . and a stochastic centralized approach has 
been pursued in fl9l . 

Apart from exploiting the reactive power capabilities of 
smart inverters, active power curtailment has been advocated 
as an ancillary service as well ll20ll . ITill . Il22l . Il23l . Droop- 
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based active power curtailment has been proposed as an effi¬ 
cient means for overvoltage prevention (20) . In ETl . an SDP- 
based relaxation has been devised for jointly commanding 
active and reactive power setpoints to inverters in multi-phase 
distribution grids. Leveraging joint reactive power compensa¬ 
tion and active power curtailment, a multi-objective OPF is 
formulated for unbalanced four-wire distribution grids in 122) . 
Local voltage control strategies are developed for customers 
enrolled in a feed-in tariff (FIT) policy 1251 . An FIT power 
supply policy compensates DG owners for the surplus of 
renewable energy they inject into the grid. Similarly structured 
programs have been successfully deployed in Europe and 
several US states l24l . 

Existing energy management schemes enforce inverter- 
related and voltage regulation-related constraints at all times. 
However, the operation of future grids could benefit from 
currently unexploited system flexibilities. Two possible options 
are the overloading tolerance of inverters and the voltage 
regulation margins. Specifically, the inverters found in DG 
units and storage units are manufactured to operate higher 
than their nameplate apparent power rating (25) . Actually, this 
feature has already been exploited in designing panels (26). 
Moreover, most voltage regulation standards, such as the 
American National Standard Institute (ANSI) C84.1 (27l . and 
the EN 50160 standard (28), define two voltage regions: 
one for normal operations, and one whose use is limited in 
frequency and duration. According to the EN 50160 standard, 
for example, nodal voltage magnitudes are required to lie in 
the latter region for 95% of any 10-minute sample (5), (28) . 

To exploit such flexibilities, this work proposes an energy 
management scheme, where voltage regulation and inverter 
capacity constraints are imposed in a stochastic rather than a 
deterministic sense. Our contribution is not on the effect of 
pricing and curtailment policies on renewable integration. It 
is rather an algorithmic framework for exploiting the afore¬ 
mentioned sources of flexibility to lower costs and improve 
renewable integration. Different from existing schemes, op¬ 
erational constraints are relaxed instantaneously, while tighter 
limits are enforced in an average sense. This is achieved using 
a stochastic dual subgradient scheme that relies on power 
flow models with different accuracy-complexity trade-offs. 
The schemes are based only on data to command set-points to 
DG units, and enjoy convergence and feasibility guarantees. 
Numerical tests using synthetic and real data on a 56-bus and 
the IEEE 123-bus grids corroborate the efficacy of the scheme. 

Paper outline. The rest of the paper is outlined as follows. 
The novel energy management problem is formulated in 
Section El An ergodic optimization approach is presented 
in Section ED while a stochastic approximation solver is 
developed in Section [Iv] The implementation of the solver 
depending on two grid models is presented in Section [VJ while 
performance advantages over the deterministic alternatives are 
supported by the numerical tests of Section [VlJ Concluding 
remarks are drawn in Section [VTlI 

Notation. Lower- (upper-) case boldface letters denote col¬ 
umn vectors (matrices), with the exception of power flow 
vectors (P,Q). Calligraphic symbols are reserved for sets, 
while R x denotes the set of all nonnegative TV-dimensional 


vectors; the symbol T stands for transposition; and 0 and 1 
denote the all-zeros and all-ones vectors, respectively. 

II. Problem Formulation 

Consider a distribution grid comprising N + 1 buses. The 
grid is modeled by a tree graph T := ( J\f 0 ,C ), where 
AT„ := {0,1,..., N} is the set of nodes (buses) and |£| = N 
denotes the cardinality of the edge (line) set C. Note that albeit 
structurally meshed, distribution grids are usually operated as 
radial. The tree is rooted at the substation bus indexed by n = 
0, and all non-root buses comprise the set M := {1,..., N}. 
Let v n be the squared voltage magnitude at bus n, and p n +jq n 
the complex power injection at bus n for all n £ J\f 0 . For 
notational brevity, nodal quantities related to non-root buses 
are stacked in column vectors v, p, and q. 

Active and reactive power injections at bus n are split 
into their generation and consumption components as p n := 
Pn ~ Pn an d Qn '■= Qn ^ Qn- For a purely load bus, the 
consumption components are oftentimes related via 

a constant power factor, whereas p® = g® = 0. A DG bus 
in addition to the nonnegative components (p°,g°), it also 
provides renewable generation p® > 0 and reactive power 
support <y®. All buses are henceforth assumed to be constant 
power buses. For buses having only a shunt capacitor, it holds 
that p® = pjj = = 0 and g® > 0. Generation and 

consumption components are stacked accordingly in vectors 
p c , p®, q c , and q®. 

The energy management controller is run centrally by the 
utility and communicates set-points to DG units. Although co- 
ordinative control of power inverters and utility-owned voltage 
regulating devices would only improve performance, it is a 
nontrivial task and is not considered here; see HE) for a dy¬ 
namic programming approach. In the envisioned scenario, the 
grid operation is divided into short control periods indexed by 
t. The duration of these periods depends on the variability of 
active and reactive power consumption and the availability of 
data predictions. Since power inverters provide a continuously- 
valued control variable that can be adjusted in milliseconds, 
transients have been reported to be negligible (8). This is 
in contrast to conventional voltage regulating equipment that 
results in switching transients. Without loss of generality, a 
30-sec control interval will be presumed hereafter. 

During time period t, the grid operator can buy or sell 
energy po,t from or to the main grid through the substation bus 
via the real-time market. The price for this energy exchange is 
7To,t, and it is assumed to be positive at all times. Apparently, 
if the real-time market operates on a 5-min basis, the price 
7To,t remains constant over 10 consecutive control periods. 
Internally in the distribution grid, customers with renewable 
generation units, such as PVs or wind micro-turbines, can 
subscribe to a so-termed FIT program; see e.g., (24), (29) . 
According to this program, the surplus of renewable energy a 
customer may have can be bought if deemed appropriate by the 
utility company at the FIT price 7r f t - Although FIT prices are 
currently adjusted on a monthly basis, time-varying prices ttf j 
are considered here for the sake of generality. Feed-in-tariff 
prices are also assumed to be positive. Energy consumption 
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from both FIT and regular customers is charged at a retail 
price 7r r j . The energy cost for customer n during period t is 
the product of 

Kr,t\Pn,t ~ Pn,t]- ~ ^u\Pn,t ~ Pn,t]+ (!) 

times the duration of the control period, where the operators 
[a]-)_ := max{a, 0} and [a]_ := max{0, — a } lf30l . Apparently, 
at most one of the two summands in ([U is nonzero per slot t. 

From the utility side, the energy cost for time slot t is 

TTO.i PO,t (P?, qf ) + T^/.t 1 T [P? - Pt ] + - TTr,11 T [P? - Pt ] - (2) 

multiplied by the slot duration, where [a] + and [a]_ for 
vectors are applied entry-wise now. Heed that the energy 
exchange with the main grid Po,t(pf ; q®) depends on the 
internal consumption and generation, and the associated power 
losses on distribution lines. Thus, the energy exchange po.t can 
be thought of as a function of the control variables (pf,q®), 
while its dependence on (p£,qj) and grid power losses is 
implicitly indicated by the subscript t. 

If the energy management scheme were to minimize the 
utility’s cost in ©, it would force the minimum possible local 
generation. To see that, consider a node n where at time t 
the demand is higher than the installed solar capacity; hence, 

-[Pn,t - Pn,t]~ = Pn,t ~ Pn,t < Then - the Utilit y EMS 
would command p 9 nt = 0 unless there is an under-voltage 
condition. Such a policy contradicts the purpose of an FIT 
program. The FIT program should curtail renewable power 
only if a customer has a surplus and the surplus cannot be 
bought due to either financial or voltage regulation reasons. 
To accommodate the FIT terms, the utility does not curtail 
renewable generation when the net injection is negative. Thus, 
the cost to be minimized by the energy management scheme 
is 7r 0 ,tPo,t(p?,qf) + 7r/,tl T [pf -p*]+ rather than that in ©. 

Operation of the energy management scheme is detailed 
next. Before control period t, the EMS collects predictions for 
prices (no,t, 7r/,t), loads (p£, q£), and the maximum renewable 
generation pf. At every period t, buses are partitioned to those 
having a renewable energy surplus comprising the set 

Aft := {n G Af : p 9 n , t > p^ t } (3) 

and to those buses belonging to the complementary set of Aft 
denoted by Aft ■ To jointly perform active power curtailment 
and reactive power management, the EMS could solve the 
ensuing problem per time interval t 


min 

p® .n? 

7ro,t po,t(Pt, q *) + nf,t i T [pf - P t}+ 

(4a) 

s.to 

o <Pn,t <P 9 u,t, V neAft 

(4b) 


Pn,t =P 9 n,t, VneFt 

(4c) 


(Pn,t) 2 + ( Qn,t ) 2 < Sn, Vn 

(4d) 


Vi < V n ,t(Pt,<lt) < v u , V n. 

(4e) 


Power injections {(p® t ,g® t )} n are constrained in the fea¬ 
sible set defined by (Bbl i-fiell. Constraints (I4bl >- (l4dt apply 
locally per bus n, whereas the voltage constraints in (Bel l 
couple power injections across the grid. Specifically, the term 
Pnt~Pnt i n is the active power curtailed for all inverters 
with renewable surplus at time t, i.e., n £ Aft- Constraint 


( l4dl > originates from the maximum apparent power capability 
(nameplate rating) s n of inverter n. Constraint (l4el > maintains 
the squared voltage magnitudes in the prescribed interval 
V := [vi, Vy\ at all nodes. Similar to the energy exchange 
Po,t(Pt,qt)> voltage magnitudes are expressed as implicit 
functions of (p®,q®), whose actual function forms depend 
on the postulated grid model and are elaborated in Section [V] 
To simplify the exposition, constraints on the apparent power 
flows on distribution lines have been ignored; such limits can 
be readily incorporated using the grid model of Section IV-BI 
It is worth mentioning that policy scenarios where the utility 
accepts any energy surplus as soon as grid constraints are 
satisfied can be captured by simply setting FIT prices 7r j <t 
to zero for all t in ©. 

Problem © guarantees that all power and voltage con¬ 
straints are satisfied at all times. Nevertheless, future distri¬ 
bution grids will afford flexibilities that can be leveraged to 
lower operational costs and better integrate renewables. Two 
possible sources of flexibility are the overloading capability of 
smart inverters and the voltage regulation ranges. Regarding 
the former, a grid-tied power inverter can exceed its apparent 
power capacity for a short period of time. Indeed, power 
electronics are empirically designed to operate at even 1.2-1.5 
times higher than their nameplate rating f25l . For the latter, 
instead of requiring the squared voltage magnitudes to lie in V 
at every t, it suffices for their time-averages to lie in V, and the 
instantaneous values to lie within a wider range. For instance, 
according to standard EN 50160, voltages are required to stay 
in V for 95% of any 10-minute sample l28l . Additionally, 
heed that problem © depends on predictions (p£, qj', p®), 
and prices {tto,u 7r /,t)- E ' s therefore optimal only if load 
demand, renewable generation, and prices are perfectly known. 
In practice though (p£, q£, p®, 7To,t, 7r/ it ) involve uncertainties 
(e.g. measurement noise, time-delay, and system variability) 
rendering the solution of © hardly optimal. 

To leverage operational flexibilities and cope with un¬ 
certainties, a stochastic rather than the deterministic energy 
management formulation of © is pursued next. To that end, 
the time-varying problem parameters {p£, q£, pf, 7To,t, 7r/,t} 
are modeled as stochastic processes 1311 . l32l . lf33l , To capture 
ensemble averages via time averages, the aforementioned 
stochastic processes are assumed stationary and ergodic, yet 
not necessarily independent across time; see |34] and [35]. 
Recall that a stochastic process is ergodic if its moments (e.g., 
the mean) can be inferred from a single realization of the 
process. The novel energy management scheme entails solving 
the following stochastic optimization problem 


J 2~ f min, E [7To,tPo,t(P?,q?) + 7r /jt l T [pf 

s.to 0 <p 9 n t <p 9 n ^ V n <E Aft 
P 9 n,t =P 9 n,t, Vnel, 

« t ) 2 + « t ) 2 <^, Vn 
Pi < Wn,t(p?,q?) < Vu, V n 
E [( Pn,t ) 2 + (<t) 2 ] <sl, v n 
Vi <E [w n ,t(pf,qf)] < v u , V n 


-P t}+] (5a) 

(5b) 

(5c) 

(5d) 

(5e) 

(5f) 

(5g) 


where the optimization variables consist of (pf,qf) for all 
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periods t, and the expectations are taken over the joint distri¬ 
bution of (pj, q£, p®, 7To,t, 7T/,t) across all periods t. Constraint 
(TSfb guarantees that the average apparent power complies with 
the nameplate inverter capacity for all buses; while constraint 
©D enforces a hard limit s n (s n > s n ) on the instantaneous 
apparent power for all n. Similarly, the averages of squared 
voltage magnitudes are maintained in V according to ( [5g| , 
whereas constraint (Bel) ensures that their instantaneous values 
lie in a region V' := [u;, v u ] with V C V'. For example, 
the ANSI C84.1 requires voltage magnitudes to lie within 
V = [0.95 2 ,1.05 2 ] per unit (p.u.) of normal operation, but 
within V' = [0.917 2 ,1.058 2 ] p.u. over short durations |27l . 

Let us compare the solution of © to the minimizers 
obtained from © at every time t. Note that constraint (l4dt 
implies constraints (l5dt and (Bft . but not the converse. Like¬ 
wise, constraints Bel l guarantees ( Bel l and ( [5g| . Therefore, 
the stochastic scheme in © constitutes a relaxation of the 
deterministic problem © solved over time t. As such, the 
minimizers of © could yield a lower average operational cost, 
i.e., .1*2 < E[./|* v t ], where the expectation is taken over time /:. 

The stochastic problem in © involves infinitely many 
variables {p^qflt. Nodal power injections at time t should 
satisfy the instantaneous constraints (I5b1 >— (f5el>. Further, the 
infinitely many variables are coupled across time via the 
objective function and the average constraints ©} and (5gi, 
hence challenging the solution of ©. A stochastic optimiza¬ 
tion approach for tackling © is pursued in the next section. 


The dual function for problem © is the minimum of the 
Lagrangian function over all primal variables. Due to the 
linearity of the expectation operator, the minimization and the 
expectation operators can be interchanged. After rearranging 
terms, the dual function is thus expressed as 

N 

9(v,i,€) :=E[St("4>?)] - (^4 “ i n v l + £ n v u ) 

n=1 

where functions £,£) are defined as 

9t(v,Z,£) ■= min {*o,tPo,t(Pt > q ?) + 7r /,t lT [P? - Pt] + 

+ J2 Un [( Pn,t f + ( 9n,t ) 2 ] 
n =1 
N 

+ H(?n-l l K,t(Pi> q t)} ( 7 ) 

n= 1 

and the feasible set fl t is given by the instantaneous constraints 
in © as 

:= {(pf,qf) satisfying ©3) - (Belli. (8) 

The dual problem is obtained by maximizing the dual 
function over the dual variables, that is 

3( l/ \£*>?*) := max g{ (9) 

",€,£>o 


III. Ergodic Energy Management 

The goal of ergodic energy management (EEM) is 
to design algorithms that sequentially observe predictions 
(Pt-qtiPt , 7To,t, TTf,t), and solve near optimally the stochastic 
problem in ©. The EEM is inspired by related ideas from re¬ 
source allocation in wireless communication networks, where 
due to propagation channel uncertainties and variabilities, one 
prefers to optimize the average rather than the instantaneous 
system behavior Il36l . l37l . The key assumption is that only 
realizations of those stochastic processes are available, while 
their joint probability density function is typically unknown. 

Since optimization variables, henceforth collectively de¬ 
noted by x := ({pj,qf} t ), are coupled via expectations, 
constraints ( BT1 ) and ( [5gj i are dualized. Let v e R+, £ £ M+, 
and £ £ E x denote the dual variables corresponding to ©}, 
and the lower and upper voltage bounds in ( |5g| ), respectively. 
All other constraints are kept explicit. Using these definitions, 
the Lagrangian function of © is readily written as 

:= E [7r 0 , t po,t(pf,qf)+7r/,tl T [pf - P?]+] 

N 

+ ^ V n {E [(Pnj) 2 + ( qn,t ) 2 ] _ S n} 

n— 1 
N 

+ L { Vl ~ E K.*(p*'> q t)]} 

n— 1 
N 

+ 5Z?«^ 1E [ 1 ’"4( p f. q ?)] -M- ( 6 ) 

n—1 


Evaluating g( iz, £,£) requires solving infinitely many prob¬ 
lems of the form in ©, and then averaging the optimal 
costs over the joint probability density function (pdf) of 
{Pt, q t.Pt , 7To,t, 7r/,t}- Even if the joint pdf were available, 
finding the expectations would be non-trivial. Hence, even 
evaluating the dual function becomes challenging. To max¬ 
imize the dual function in a feasible manner, a stochastic 
optimization solver is proposed next. 


IV. Stochastic Approximation Solver 

The problem at hand is tackled using a stochastic dual 
subgradient method Il36ll . l37l . l38l . To maximize g(u. £, £), 
the Lagrange multipliers are updated using the projected 
subgradient iterations for some step size /i > 0, as 


v t := [i/ t _i 

t := 


■Mj.t 


J + 


: = €t -1 + l jS 7. 




(10a) 

(10b) 

(10c) 


where the vector S t := [S~l t dj t 6j 


is a subgradient of 
gt(y, evaluated at the previous iterate (vt-i, £ t _ 1 i £t-i)- 
The entries of the subgradient vector, denoted by can 

be found as 



■■=(P 9 n,t) 2 + K t ) 2 -S 2 n 

(11a) 


■= Vi - tvt(pf,qf) 

(lib) 


:= u„,t(p?,q?) - v u 

(He) 
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TABLE I 

Ergodic Energy Management Algorithm 


1: Input operational limits {s„,s„}„ 6 u, {vi,v u ), (Hi,v u ), 
and step size p > 0. 

2: Dual variables zz 0 , £ and £ 0 are initialized to zero. 

3: For t = 1, 2,... do 

4: Acquire predictions (pt, q?, p?, 7To,t, 7r/,t). 

5: Find primal variables (p®,qf) as the minimizers of 
b y solving C© or <2Qj. 

6: Update dual variables (i' t . , £, t ) using Go). 

7: Communicate setpoints (p®,q®) to DGs. 

8: End for 


for all n, where (p®,q®) are the minimizers of the problem 
in © for g t (ist-i,£ t ,Note that the Lagrange multi¬ 
pliers are updated at every control interval. 

Table |H summarizes the EEM algorithm. Operational limits 
as well as the step size are set in Step 1, and Lagrange 
multipliers are initialized to zero in Step 2. The EEM then 
iterates between four steps. In Step 4, the utility collects 
predictions for the random variables involved. In the absence 
of more elaborate options, the most recently observed or 
metered values can be used as predictions for the upcoming 
control period of interest. Step 5 finds the optimal primal 
variables by solving © evaluated at the current value of the 
Lagrange multipliers. Step 6 updates the Lagrange multipliers 
via the dual subgradient rule of (ITOl) . The calculated setpoints 
(p®,q®) are finally communicated to the DGs, and applied 
to the grid in Step 7. It is worth stressing that the proposed 
EEM scheme does not require any distributional knowledge on 
the input data (p£, q£ , p®, 7To,t, Moreover, although the 
focus is on utility cost minimization, other energy management 
tasks such as voltage regulation and conservation voltage 
reduction, could be cast under this framework. 

As far as convergence is concerned, note first that 
all primal and dual iterates depend on the realizations 
(Pt.q^Pt jTTo,t)7r/,t)> and are thus random. For that reason, 
convergence claims are in probability. Using the definition 
of S t , it is easy to show that there exists a finite H such 
that E[||<5 t ||||iz f _ 1 ,| ( _ i ,f t _ 1 ] < H 2 for all t, i.e„ the 
subgradient 6 t is bounded at all times. In particular, it holds 
that [S Uit \n < s 2 n , while and are both upper 

bounded by (v u — vf) 2 . Thus, the bound H can be selected as 


1 v * i 

vi < lim - r (p?,q?) < v u 

£—>•00 t z ' 


(13b) 


T = 1 


Furthermore, the incurred operational costs satisfy 

1 


pro,- 


•P0,r(p®,q®)+7T/, 


-l T [p?- 


PrJ + 


— J< 


2 — 


pH 2 


t 

lim 4 
t—foo t J 
r=1 

almost surely for H as in ©. 

The proof of Proposition |T] can be found in lfT7l . Proposi- 
tion|T|asserts that the ensembles of primal sequences { p® . q®} t 
are feasible almost surely, meaning that constraints (f5fb and 
(5g l are satisfied almost surely. Moreover, the ergodic limit of 
the objective is at most pH 2 /2 away from the optimal J% [cf. 
©]. The aforementioned claims hold even if the stochastic 
processes involved are correlated across time 1371 . Although 
stochastic processes have been assumed to be ergodic for the 
theoretical claims to hold, the numerical tests in Section [VI] 
using real data show the efficacy of the scheme even with 
non-ergodic data. 

The EEM problem in © and its stochastic approximation 
solver of Table Q] involve the power losses po,t(p?)Qt) an d 
the squared voltage magnitudes {v njt (pf, qf)}„. So far, both 
quantities have been expressed as functions of the control 
variables (p®, q®). In that respect, the EEM scheme constitutes 
a general framework where different power system models can 
be assumed. To implement Step 5 in the algorithm of Table © 
the actual forms of po,t(Pt , q?) and {u n , t (pf, q®)}„ need to 
be specified. As a turnkey application of EEM, the ensuing 
section focuses on radial single-phase distribution grids using 
two power flow models with different accuracy-complexity 
trade-offs. 


V. Grid Modeling and Algorithms 

This section specifies functions Po,t(Pt. Q?) an d 
{vn : t(Pt ■ q?)}n using an exact full AC grid model and 
its linear approximation. Both cases are then integrated into 
the EEM algorithm. Selecting between the two models relies 
on the computational capabilities that can be afforded. The AC 
model-based EEM can be formulated as an SOCP, whereas 
the linear model yields a linearly constrained quadratic 
program. Therefore, the latter option offers an approximate 
yet computationally less demanding alternative to the former. 

A. Branch Flow Model-based EEM 


N 

# :=EK + 2 (^-^) 2 ]- < 12 ) 

n— 1 

Adopting l37l Theorem 1], the following result characterizes 
the almost sure feasibility and optimality of the EEM algo¬ 
rithm. 

Proposition 1 (1771). For the sequences {p®,q®}t generated 
by the algorithm in Table 0 the next hold with probability 1 
for all n G Jf 

Km \ J2^P 9 n,r) 2 + (€,r) 2 } < 4 


Due to the radial structure of distribution grids, every non¬ 
root bus n G Af has a unique parent bus, which will be denoted 
by a n . The directed edge ( a n ,n ) G C corresponding to the 
distribution line feeding bus n will be indexed by re, see Fig. Q] 
Without loss of generality, buses can be indexed such that 
a n < n for all n gAI. 

Let z n = r n + jx n be the line impedance of line n, and 
£ n j the squared current magnitude on line n at time t. If 
Pn,t +jQn,t is the complex power flow on line n seen at the 
sending end bus a n at time t, the grid can be described by 
the branch-flow model |f39l 

Pn,t = ^ ^ Pk,t ( Pn,t T n f n ,i) 

fcGCn 


(13a) 


(14a) 
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Fig. 1. Bus n is connected to its unique parent a n via line n. 


Qn,t — ^ ^ Qk,t ( Qn,t 3Cn@-n,t) 

keC n 

Vn,t — ^(y.n,t ‘2j{y'nF > n,t %nQn,t') H - {j*n 
c _ P n,t + Qn,t 

n.t — 

Voi n ,t 


(14b) 

(14c) 

(14d) 


for all n £ AT, where C„ is the set of the children nodes 
of bus n. The power injections at the substation bus are 

Po,t = J2k&c 0 Pk ’ u 9o,t = ^fceCo Qm> and its squared 
voltage magnitude no, t is controlled at a desirable value. 
Similar to (p*,qt). the vectors r, P t , Q t , v t , and £ t , collect 
the entries of r n , P n , t , Q n ,t, v n ,t, and £ n>t , accordingly. 

Equations (1 1 4ab - (l 1 4cb are linear with respect to the system 
variables (p t , q t , P f , Q t , v t , £ t ). The equations in (1 1 4db are 
quadratic yielding a non-convex feasible set. Nonetheless, 
the latter equations have been recently relaxed to convex 
inequalities described by the hyperbolic constraints in 


P n,t Qn.t — v a n ,t£-n,ti V 


(15) 


which can be equivalently expressed as the convex second- 
order cone constraints 


2 P n ,t 
2Qn,t 

£n,t 


— VfXrnt 


4“ £-n.ti ^ ti. 


(16) 


Equations (1 1 4ab - (l 1 4cb and (fl6l > represent now a convex feasi¬ 
ble set. Recent works suggest using this relaxed feasible set 
to perform several grid optimization tasks. Under different 
conditions, the relaxation has been shown to be exact, i.e., the 
obtained minimizer attains ( 1 1 6[ i with equality; see EL and ref¬ 
erences therein. Henceforth, all SOCP relaxations are assumed 
exact, which will be numerically verified in Section [VT] 
Based on d!4aM14ct and (fl6l >, the active power injection 
at the substation bus po,t(Pt ? Q?) can be expressed as 


N 


N 


PoAPt, qf) = ~ Pn,t) + 


= l T (p?-pf)+r T ^ t 


where the second summand represents the total power losses 
on distribution lines. Hence, under the aforementioned relaxed 
grid model, the primal update (Step 5 in Table Q} entails solving 
the optimization problem 

min 7r 0 , t l T (Pt - pf) + 7r 0 , t r T £ t + 7r/, t l T [pf - pf] + 
p“.q; At, 
p t -Qt > v t 

+ "n,t -1 [{Pn,t) 2 + KA 2 ) 

n= 1 


N 

+ 1 ~ t_ 1 ) V n,t (17) 

n= 1 

s.to (l5bl) — (f5el). (I14ab — (U4cb . (fl6l>. 

In addition to the original variables (pf, qf), the primal up¬ 
date now involves the variables (P t , Q f , v t , £ t ) too. Problem 
(fT7l can be reformulated to an SOCP. All instances of (fl7l> 
solved in Section [Vi] were exact. Nevertheless, solving (fTTb 
could be computationally demanding for large-scale distribu¬ 
tion grids. This motivates our next instantiation of the EEM 
algorithm under an approximate grid model. 


B. Linear Distribution Flow-Based EEM 

The linear distribution flow model can be derived as follows. 
Because the line parameters {r n , x n } ne j\f have relatively 
small entries, the last summands in (1 1 4ab -( IT4c1) can be ignored 
yielding the linear equations for all n £ M [|39l 


Pn,t — ^ ^ Pk,t Pn,t 

(18a) 

k£C n 

Qn,t — ^ ^ Qk,t Qn,t 

(18b) 

kec n 

Vn,t — 'V C t n ,t e ^{j'nPn,t H"~ %nQn,t) • 

(18c) 


In this way, squared voltage magnitudes are now approximated 
as linear functions of (pt, qt). Assuming squared voltage mag¬ 
nitudes to be close to 1 p.u., squared line current magnitudes 
are approximated as ll39l 


P 2 

n,t 


■Qn,t 




Qr, 


(19) 


" OL n ,t 


Therefore, the active power injection at the substation bus can 
be thus approximated by 


N 

PoAPt, qf) = ! T (Pt - pf) + Vn ( P lt + Qn.t) • 

n—1 

Building on the approximate model of (fl8l)—(IT9t . the primal 
update of the EEM algorithm (Step 5 of Table Q} for period t 
entails solving the problem 


N 

min 7T 0 , t l T (Pt - pf) + 710,4 r A P n,t + Ql,t) 

Pt ’ q * ’ nA 

N 

+ 7T/,tl T [pf - Pf]+ + “ in,t-l> Vn ^ 

n= 1 
N 

+ 1 [(Pn ,t) 2 + AnA 2 ] ( 20 ) 

n= 1 

s.to (f5bl> — (f5el). (1 1 8ab — (1 1 8cb . 

From (fl8all-(fT8bl). the line flow variables (P t ,Q t ) can be 
substituted as linear functions of (pt, qt). Hence, problem ( 120| ) 
can be solved as a linearly constrained quadratic program. 
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TABLE II 

Line data for the 56-bus distribution feeder fTTli 


From 

Bus 

To 

Bus 

n 

P] 

Xi 

P] 

From 

Bus 

To 

Bus 

n 

P] 

Xi 

P] 

i 

2 

0.160 

0.388 

28 

29 

0.395 

0.369 

2 

3 

0.824 

0.315 

29 

30 

0.248 

0.232 

2 

4 

0.144 

0.349 

30 

31 

0.279 

0.260 

4 

5 

1.026 

0.421 

32 

33 

0.263 

0.073 

4 

6 

0.741 

0.466 

32 

34 

0.071 

0.171 

4 

7 

0.528 

0.468 

34 

35 

0.625 

0.273 

4 

20 

0.138 

0.334 

34 

36 

0.510 

0.209 

7 

8 

0.358 

0.314 

34 

38 

1.062 

0.406 

8 

9 

2.032 

0.798 

34 

41 

0.115 

0.278 

8 

10 

0.502 

0.441 

36 

37 

2.018 

0.819 

10 

11 

0.372 

0.327 

38 

39 

0.610 

0.238 

11 

12 

1.431 

0.999 

39 

40 

2.349 

0.964 

11 

13 

0.429 

0.377 

41 

42 

0.159 

0.384 

13 

14 

0.671 

0.257 

41 

47 

0.157 

0.379 

13 

15 

0.457 

0.401 

42 

43 

0.934 

0.383 

15 

16 

1.008 

0.385 

42 

44 

0.506 

0.163 

15 

17 

0.153 

0.134 

42 

45 

0.095 

0.195 

17 

18 

0.971 

0.722 

42 

46 

01.915 

0.769 

18 

19 

1.885 

0.721 

47 

48 

1.641 

0.670 

20 

21 

0.251 

0.096 

47 

49 

0.081 

0.196 

20 

23 

0.225 

0.542 

49 

50 

1.727 

0.709 

21 

22 

1.818 

0.695 

49 

51 

0.112 

0.270 

23 

24 

0.127 

0.542 

51 

52 

0.674 

0.275 

23 

25 

0.284 

0.687 

51 

53 

0.070 

0.170 

25 

26 

0.171 

0.414 

53 

54 

2.041 

0.780 

26 

27 

0.414 

0.386 

53 

55 

0.813 

0.334 

26 

32 

0.205 

0.495 

53 

56 

0.141 

0.340 

27 

28 

0.210 

0.196 






TABLE III 

BUS DATA FOR THE 56-BUS DISTRIBUTION FEEDER HU 


Load Data 

Load Data 

Bus 

Peak 

Bus 

Peak 

Bus 

Peak 

No. 

[MVA] 

No. 

[MVA] 

No. 

[MVA] 

3 

30 

25 

0.20 

43 

1.34 

5 

0.67 

27 

0.13 

44 

0.13 

6 

0.45 

28 

0.13 

46 

0.67 

7 

0.89 

29 

0.07 

47 

0.13 

8 

0.07 

31 

0.13 

48 

0.45 

9 

0.67 

32 

0.27 

50 

0.20 

10 

0.45 

33 

0.20 

Shunt Capacitors 

11 

2.23 

34 

0.27 

Bus 

Nameplate Capacity 

12 

0.45 

35 

0.45 

No. 

[Mvar] 

14 

0.20 

36 

1.34 

19 

0.6 

16 

0.13 

37 

0.13 

21 

0.6 

17 

0.13 

38 

0.67 

30 

0.6 

18 

0.20 

39 

0.13 

53 

0.6 

19 

0.45 

40 

0.45 

Base Information 

33 

2.23 

41 

0.20 


Phase = 12kV 

24 

0.45 

42 

0.45 


Sbase = 1MVA 


VI. Numerical Tests 

The novel schemes were numerically tested on a 56-bus 
distribution grid from Southern California Edison (SCE) and 
the IEEE 123-bus feeder | jTT| , fl40j . Line and bus data for the 
SCE grid are listed in Tables Hand m accordingly, while 
a power factor of 0.8 is assumed for all loads; see [HU for 
details. The capacity of the PVs installed on buses 19 and 45 
was set to 6 MW. At each 30-sec control period, the EEM 
controller collects power demands from load buses and solar 
generation predictions from PV units. Subsequently, active and 
reactive power injections by PV inverters are determined by: i) 
solving the deterministic energy management (DEM) scheme 
in ©; and ii) the novel EEM algorithm of Table U that is 
initialized to zero. 


—«— EEM (AC, Cost=-$589) 
——- EEM (AC, average) 

—*— DEM (AC, Cost=-$565) 



-700 - 1 - 1 - 1 - 1 - 1 - 

0 10 20 30 40 50 60 

Time [min] 

Fig. 2. Energy management cost using the AC model on synthetic data 
{H = 0.1 for EEM). 



Time [min] 

Fig. 3. Energy management cost averaged over 20 independent realizations 
using the AC model. 


The margins for squared voltage magnitudes are set 
as [vi, v u ] = [0.9604,1.0404] p.u. and \v t , v u ] = 

[0.9409,1.0609] p.u., with nominal voltage vq = 1 p.u. The 
apparent power capability for smart inverters is set to 1.3 times 
the nameplate capacity of the associated PV. Performance is 
evaluated in terms of the energy management cost and the 
instantaneous counterpart of the cost in ©. All algorithms 
were implemented using MATLAB and CVX on an Intel 
CPU @ 3.4 GHz (32 GB RAM) computer BT1 . Every run 
for the full AC and the linear approximation model-based 
algorithms on the 56-bus grid was completed in 1.5 and 1.3 
seconds, respectively. The related times for the IEEE 123- 
bus feeder increased to 4.5 and 3 seconds, respectively. It is 
worth mentioning that all SOCP relaxations encountered in 
the ensuing experiments were feasible and exact. 

To verify the almost sure optimality, the first experiment on 
the 56-bus grid simulates synthetic load consumption and solar 
generation as p) = p c +e£ and pj / = p 7 respectively. The 
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nominal values p c and p 7 are set to 40% of the peak demand 
values and 80% of the maximum PV generation, accordingly. 
Vectors e£ and e 9 t capture fluctuations modeled as independent 
zero-mean Gaussian vectors with standard deviations equal to 
5% of the corresponding nominal values. Given that current 
FIT prices change on a monthly basis and are oftentimes half 
of consumption prices J23J, prices were set to 7To,t = 30(j/kWh 
and 7 Vf t t = 1 of:/kWh for all t. 

Using the branch flow model. Fig. [2] depicts the energy 
management cost for the deterministic and the ergodic energy 
management schemes over a single realization of 120 control 
periods. The step size for the ergodic scheme is set to /x = 0.1, 
while the time-average energy management cost per time slot t 
is defined as \ £‘=1 [7r 0jT p 0 ,r (p?, q?) + 7 T/, r l T [p®-p^]+] . 
The actual total operational cost over an hour is —$565 for 
the DEM and —$589 for the EEM scheme. 

The second test studies the effect of the step size fi on 
the convergence of EEM. The AC-based EEM scheme was 
simulated for fi G {0.1, 0.2, 0.3} along with the DEM scheme. 
Twenty Monte Carlo system realizations over 60 control 
periods were averaged for each step size value, while the 
corresponding average energy management cost is plotted in 
Fig. [3 The curves demonstrate that larger step sizes incur 
higher energy management costs, an observation that agrees 
with the optimality gap of Proposition [Q 

To test the proposed schemes in real-world conditions, the 
ensuing two experiments entail real data from the Smart* 
project (21. Consumption data involved the electricity usage at 
minute-level samples from 443 homes on April 2, 2011; and 
the power output of 3 residential PVs collected at 5-second 
intervals over August 12, 2011. Data were preprocessed as 
follows. Consumption data were first linearly interpolated to 
yield 30-sec loads, and then averaged over every 10 homes to 
better resemble bus loads. Daily load curves were subsequently 
normalized to a maximum value of one, and mapped to differ¬ 
ent buses DU- Normalized daily load curves were multiplied 
with the nominal load value per bus. Concerning PVs, 5-sec 
data were aggregated to 30-sec data. Daily generation data 
were likewise scaled to match rated capacities. 

A single system realization was simulated over the 600 30- 
sec control periods during 9:30am-2:30pm for both the AC- 
and the approximate model-based schemes. Figure [4] presents 
the cost for fi = 0.25. Using either the AC or the linear 
approximation model, the novel EEM scheme achieves a lower 
cost than the DEM one. 

Fig. E| depicts the evolution of the squared voltage mag¬ 
nitude for bus 45, and the evolution of the dual variable 
£ 45 ,t for the tight voltage margin constraint in (5gl. The 
voltage magnitude for the deterministic scheme remains in 
the tight region [vi, v u \ = [0.9604,1.0404] throughout the 
operation horizon. The voltage magnitude obtained from the 
stochastic scheme lies occasionally beyond the voltage margin 
v u = 1.0404. Nonetheless, over-voltage effects have short 
duration. At around 10:25 am, when the voltage magnitude 
violates the tight voltage constraint for the first time, the dual 
variable becomes positive and starts increasing. As long as 
the voltage magnitude fluctuates above the tight margin, the 
dual variable keeps increasing. After roughly 12:20 pm, the 


-EEM (AC, Cost=-$2797) 

- EEM (AC, average) 

-DEM (AC, Cost=-$2680) 

-DEM (AC, average) 

EEM (LDF, average) 

- DEM (LDF, average) 



11:30 12:30 

Time [hours of the day] 


Fig. 4. Energy management cost using the AC and the linear distribution 
flow (LDF) models on real load and solar generation data 0. 



Time [hours of the day] 


Fig. 5. Voltage magnitude for bus 45 and the associated dual variable ^ 45 ,t 
using the AC model-based schemes. 


voltage magnitude drops and remains consistently below the 
upper margin, while the dual variable decreases and eventually 
becomes zero for the rest of the day. 

To get a grid-level view of voltage regulation and active 
power curtailment, the top panel of Fig. [6] shows the grid- 
averaged voltage magnitude obtained via the DEM and EEM 
schemes, as well as without any control. Under no control, 
voltage magnitudes consistently exceed regulation margins. 
Moreover, the EEM scheme yields slightly higher voltage 
profile than DEM in exchange for lower operational cost. 
Similarly, the bottom panel of Fig. [6] shows the grid-wise 
solar generation curtailment incurred by DEM and EEM. 
Apparently, the DEM scheme curtails more active power than 
the EEM scheme. 

The last test involved the IEEE 123-bus feeder shown in 
Fig- 0 ED. The original multiphase system was heuristic ally 
modified to a single-phase one as described in CD: Loads 
were split uniformly over all phases. Line self-impedances 
























9 


v u = 1.03 

\ 





v A-/. 

* V 

fy f X l" 

; 'V , 4 

Wft,. 

/ 

v u = 1.02 



1 W 









- EEIV 

(AC) 

(AC) 

ontrol (AC) 


vi = 0.98 


- No c 



Vi = 0.97 

A 






9:30 10:30 11:30 12:30 13:30 14:30 



1.4 



9:30 10:30 11:30 12:30 13:30 14:30 


Time [hours of the day] 


Fig. 6. Top: Grid-averaged voltage magnitude using the AC model-based 
schemes. Bottom: Total active power curtailment over 9:30am-2:30pm. 


Fig. 7. Schematic diagram of the IEEE 123-bus test system with a PV (40j. 



Time [hours of the day] 


Fig. 8. Energy management cost evaluated on the IEEE 123-bus test system 
using the linearized model on real load and solar generation data. 


were averaged over phases, while mutual impedances were 
neglected. Closed switches were modeled as short circuits and 
open switches were ignored. Distributed loads were replaced 
by two identical spot loads at the two line ends. Transformers 
were modeled as lines with given impedances, and tap ratios 
for all voltage regulators were fixed to 1.08. A single PV 
with nameplate rating of 1.2 MW is installed at bus 114, 
which corresponds to PV penetration of about 100%. With 
vo = 1 p.u., voltage regulation bounds are chosen as [vi, u u ] = 
[0.9801,1.0201] p.u. and [u z , v u \ = [0.9409,1.0609] p.u., 
while inverters can be overloaded by 110% their nameplate rat¬ 
ing. The linearized model was adopted, and real data for solar 
generation and home loads were utilized. Each time period t, 
the prices were set to noj = 30<:/kWh and 7r/ jt = 15<:/kWh. 
Fig. [8] presents the cost over 600 30-sec control slots during 
9:30 am - 2:30 pm for the two schemes. The step size was set 
to // = 0.001. The total operational cost over the simulation 
period amounts to —$1,146 for EEM and —$991 for DEM, 
thus demonstrating the superiority of the ergodic approach in 
the IEEE 123-bus feeder. 


VII. Concluding Remarks 

This paper introduced an EEM framework. Smart inverters 
are engaged in active power curtailment and reactive power 
support in a stochastic sense. A stochastic dual subgradient 
scheme enforces tighter operational margins at all times, yet 
letting system characteristics deviate over short time intervals. 
The developed algorithms are guaranteed to converge to the 
optimal operational point, while the feasibility is satisfied 
almost surely. Numerical tests using a full AC grid model 
and its linear approximation on a 56-bus grid and the IEEE 
123-bus feeder demonstrated the viability of the approach. In 
particular, the grid was operated within the regulated margins 
at all times, while local variables could fluctuate over looser 
ranges during extreme conditions. The suggested flexible grid 
operation brings up several interesting questions. Enforcing 
probabilistic rather than average constraints is worth investi¬ 
gating. Decentralized and localized implementations are timely 
and pertinent. Integrating utility-owned voltage regulating 
equipment to develop coordinative control schemes constitutes 
an interesting and challenging future research direction. 
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